clc,clear
GR=zeros(1,22100);
for i=1:22100;
GR(1,i) = GR_x(i);
end
figure()
plot(GR)
title('GR-x')

% CP=500;
% Va=3;  %风速
% mT=90; %总质量
% alpha=1/2*1.2234*(0.2565+0.0044);
% beta=(0.0032*mT*9.80665);
% v=11*exp(0);
% p=alpha*v*(v+Va)^2+beta*v*cos(atan(GR))+v*mT*9.80665*sin(atan(GR));

% for i=1:22100
% v = find_v_Ptot(GR(1,i),500);
% end
% figure()
% plot(p)

% Va=3;  %风速
% mT=85; %总质量
% alpha=1/2*1.2234*(0.2565+0.0044);
% beta=(0.0032*mT*9.80665);
%%
CP = 400;
%% var=0 时所耗时间
V=zeros(1,22100);
for i=1:22100
A=Solve3Polynomial(0.1596,6*0.1596,(2.8243*cos(atan(GR(1,i)))+85*sin(atan(GR(1,i)))),9*0.1596-CP);
V(1,i)=A(find(imag(A)==0));
end
figure()
plot(V)
title('V-x')
t=1900./V(1,1800) + 800./V(1,2600) + 2200./V(1,4800) + 5400./V(1,10200) + ...
    5200./V(15400) + 2000./V(17400) + 4600./V(22000);

t=t./60.*2

%% 
for n=0.01:0.01:0.1
var = n;
P=ones(1,22100);
P(1,1:4899)=CP*(1-var);
P(1,4900:10299)=CP*(1+var);
P(1,10300:15499)=CP*(1-var);
P(1,15500:17499)=CP*(1+var);
P(1,17500:22100)=CP;
figure()
plot(P)
title('P-x')
for i=1:22100
A=Solve3Polynomial(0.1596,6*0.1596,(2.8243*cos(atan(GR(1,i)))+85*sin(atan(GR(1,i)))),9*0.1596-P(1,i));
V(1,i)=A(find(imag(A)==0));
end
figure()
plot(V)
title('V-x')
t=1900./V(1,1800) + 800./V(1,2600) + 2200./V(1,4800) + 5400./V(1,10200) + ...
    5200./V(15400) + 2000./V(17400) + 4600./V(22000);
t=t./60.*2
end
